In this paper, we introduce a novel strategy for efficient global robust optimization leveraging the two-stage metamodeling approach to overcome an optimization under uncertainty, which involves expensive function evaluations. In the first stage, Gaussian Process Regression (GPR) is employed to represent high-fidelity data obtained from initial designs. An adaptive infill sampling technique based on expected improvement is then used to enhance the accuracy of the global model. The output distribution resulting from the propagation of uncertainty variables is estimated using a stochastic simulator applied to the GPR model. In the second stage, Multi-Layer Perceptron (MLP) neural networks are utilized to quantify uncertain moments. To achieve an optimal robust solution, we propose a multi-objective design optimization approach using a non-dominated sorting genetic algorithm. This approach aims to balance performance mean and variation criteria for decision-making. The performance of the proposed algorithm is evaluated through a test case, demonstrating significant improvements in global model accuracy and convergence towards a reliable robust solution. We further validate this approach through its application in acoustic noise reduction via shape optimization, where consideration of uncertainty is critical.